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INTRODUCTION 



The Frobenius-Perron (FP) operator plays a central 
role in the statistical analysis of chaotic dynamical sys- 
tems by ruling the time evolution of distribution func- 
tions (probability densities) in phase space. In flows, it 
leads to a continuity differential equation known as the 
Liouville equation. The FP propagator admits a Hilbert 
space representation as a unitary operator, whose spec- 
trum must therefore belong to the unit circle. The struc- 
ture of the corresponding spectral decomposition and the 
nature of the phase space dynamics are intimately con- 
nected. Besides, this same structure determines the be- 
havior of time correlation functions and their frequency 
spectral densities. For instance, the FP operator for the 
motion on an n-torus in an integrable Hamiltonian flow 
has a purely discrete spectrum with eigenvalues given by 
e* kot , where k is any n- vector of integer numbers and 
ft is the n- vector of the torus fundamental frequencies; 
time correlation functions are in this case quasiperiodic 
and the corresponding spectral densities present delta- 
function singularities at certain frequency values from the 
set to = k$7. On the other hand, in a mixing dynamical 
system, if we exclude the eigenvalue 1 (which is simply 
degenerate and corresponds to the invariant measure), 
the spectrum is continuous; time-correlation functions 
must therefore decay. In many discrete maps and con- 
tinuous flows with chaotic dynamics this decay is expo- 
nential; however, even in this case, the correlation func- 
tions may present strong oscillatory modulations; these 
lead to bumps in the corresponding frequency spectral 
densities which have been interpreted as resonances, i.e. 
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poles in those functions after their analytical continua- 
tion into the complex frequency plane. From the original 
theoretical work carried out on this problem by Pollicott 
jl| and Ruelle Q for a class (Axiom-A) of chaotic dy- 
namical systems these singularities are generally known 
as Pollicott-Ruelle (PR) resonances. If, as happen in 
this class of systems, the location of the singularities in 
the complex frequency plane is an intrinsic property of 
the dynamical system, i.e. independent of the observable 
monitored, then these resonance poles may be considered 
to be in correspondence with generalized eigenvalues of 
the FP operator ||. This interpretation requires exten- 
sions of this operator which hold for either positive or 
negative times, and which make U a non unitary prop- 
agator whose spectral decomposition may be written in 
terms of the biorthonormal basis set of its left and right 
eigenstates. A physical way of defining the generalized 
eigenvalue problem is as an usual Hilbert space spectral 
analysis of the coarse-grained dynamics in the limit of 
zero coarse graining. This particular way establishes a 
correspondence between the coarse-grained Liouvillian 
dynamics in a chaotic flow and the classical diffusion 
equation for disordered systems j|. Such a correspon- 
dence is the origin of a new approach to carry out the 
statistical analysis of the eigenvalue spectrum of a quan- 
tum system whose classical limit is chaotic, which is an 
extension of the field theoretical methods used for quan- 
tum disordered systems ||, In this way the role played 
by the FP operator and its resonances is important not 
only in the statistical analysis of classical chaotic systems 
but also in that of their quantum counterparts ||. 

The determination of the leading PR resonances is a 
necessary step in those statistical analyses. In general, in 
order to extract them one has to resort to numerical ap- 
proximate schemes. However, the numerical approaches 
followed so far are few and not always well founded both 
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mathematically and physically. The theoretically most 
complete treatments are also the most involved numer- 
ically. The methods based on the diagonalization of a 
coarse-grained FP operator in the Hilbert space of square 
integrable functions belong to this class. In the limit of 
zero coarse-graining and a complete basis set one would 
obtain the resonances from the eigenvalues. In order to 
obtain this limit accurately, repeated diagonalization of 
a very large matrix (e.g. 8100x8100 for the simple stan- 
dard map 0) at different values of the coarse- graining 
parameter are required; in some cases this parameter is 
already set by the matrix dimension ||. A different ap- 
proach also in this class, uses a periodic orbit represen- 
tation of the Fredholm determinant of the FP operator, 
which can be expressed in terms of the traces of this oper- 
ator at different times; the location of the resonances are 
then obtained from the zeros of the corresponding Zeta 
function ||[ 0, [l(J. The scheme requires the knowledge 
of a large number of periodic orbits, which can become a 
very difficult task in many real systems. 

We already know that the decay behavior of time corre- 
lation function between two observables and the analyt- 
icity properties of the corresponding frequency spectral 
densities are intimately related to the generalized eigen- 
values of the FP operator. This is the origin of a category 
of methods in which one actually carries out an analyti- 
cal continuation of the spectral density into the complex 
frequency plane in order to find the poles associated to 
the resonances contributing simultaneously to the two 
observables chosen. Only a few examples that take this 
approach can be found in the literature. It was first fol- 
lowed by Isola for the Henon map [jll| and later by Baladi 
et al. (1^1 for intermittent systems. As will be seen in 
this paper, the variational method proposed by Blum and 
Agam also belongs to this class. This scheme, unlike 
the class of methods discussed in the previous paragraph, 
requires very small computational efforts. For instance, 
Isola Jllj] finds the leading resonances of the Henon map 
from a [L, M] Pade approximant to the spectral density 
with 16 < L + M < 26; similarly, Blum and Agam j| 
obtain the leading resonances of the cat and standard 
maps as the eigenvalues of a 4 x 4 matrix. However, the 
mathematical and physical foundations of this approach 
are not as deep as in the previous methods. This is seem- 
ingly due to the lack of a general theoretical framework 
from which one can derive not only these approximate 
numerical schemes, but also error and convergence crite- 
ria. Finding such a framework is the main goal of this 
work. It was achieved by relating these methods to the 
memory function techniques used in the general theory of 
relaxation or to the filter diagonalization approach 
followed in the location of quantum resonances |l4|, 
A new class of numerical methods to determine PR res- 
onances comes out from this theoretical framework. 

The paper is organized in the following way. In Sec- 
tion II we will present the results of the memory func- 
tion method which are more relevant to our study. These 
techniques are usually applied to Hamiltonian, Liouvil- 



lian, and general relaxation operators |U|. In determin- 
istic dynamical systems, this would imply dealing with 
the Liouville operator. However in this case the memory 
function scheme does not perform the analytic continua- 
tion required in the determination of the PR resonances 
unless a coarse grained version of that operator is used. 
If instead one chooses to analyze the FP operator this 
coarse graining is not needed; but in this case we will 
require a particular implementation of the memory func- 
tion methods to deal with propagators. This is all per- 
formed in Section II. We will show, for instance, how 
Pade approximants appear naturally in this scheme. Fil- 
ter diagonalization [ |l4| |l5[ is another approach which is 
specifically adequate to our problem. This is a particu- 
lar formulation of the harmoni c in version problem as an 
eigenvalue equation. In section III we give an account of 
this method and show its total equivalence with the mem- 
ory function scheme. The results of these two sections set 
up a theoretical framework for a new class of numerical 
schemes to determine PR resonances in chaotic systems. 
Besides, we prove in Section IV that all these procedures 
are equivalent linear formulations of the non-linear nu- 
merical problem known as interpolation by exponentials 
fljjl . This connection provides new tools to perform a 
better analysis of issues such as convergence and numeri- 
cal stability, which are relevant to establish the reliability 
of the numerical approach presented here to locate PR 
resonances. From such analysis improved schemes may 
be designed, as the one that we propose in this section 
based on a least squares fit method. The application of 
these methods is illustrated in simple dynamical systems 
in Section [v| Finally, a summary is presented in Section 

vil- 



li. THE MEMORY FUNCTION TECHNIQUE 

As mentioned in the introduction, the direct applica- 
tion of this technique to deterministic dynamical sys- 
tems would imply dealing with a coarse grained Liou- 
ville operator. If instead one chooses to analyze the 
Frobenius-Perron operator the coarse graining is not gen- 
erally needed. In this section we will proceed to perform 
the required particular implementation of the memory 
function methods to deal with this propagator. 

First of all, we have to make the following choice for 
the scalar product between observables, 
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where [i is the natural invariant measure for the dynam- 
ics. The time evolution of the system is ruled by its 
Frobenius-Perron operator U. In the case of maps this 
corresponds to one iteration step; for flows, U may corre- 
spond either to the Poincare map or to a given finite-time 
step. We then have 



|/(n + l))=f/|/(n)). 



(2) 
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Making use of this notation the autocorrelation function 
G(n) for an observable / reads 



C(n) = </(0)|/(n)) = </(0)|tP|/(0)). 
We define a resolvent G(u>) of U by writing 



(3) 



G(o;)|/(0)} = £e^|/(n)) 

n=0 

oo 

= ^(e-t/)"|/(0)), 



n=0 



Therefore, 



G( W ) 



1 - e iu) U 



Imu > 0. (4) 



(5) 



In the memory function formalism [|13| one chooses a 
particular state |/) = |/(0)) and considers the diagonal 
matrix element G//(w) = (/ | G{oj) | /). This is an 
analytic function in the upper half plane of the complex 
frequency co. Its singularities, generally in the form of sin- 
gle poles at co = C0i, can appear in the lower half plane. 
The eigenvalues of U in term of these poles are therefore 
e~ luJi . Recursive projective techniques, first developed by 
Zwanzig [|l7j and Mori Jl8[ , are then implemented to ex- 
press Gffjco) as a continued fraction. The process starts 
out by defining two complementary generalized projec- 
tion operators Pq = U^Hfel and Qo = 1 — Po, with the 

(/o|/o> 

identifications |/o) = |/) and (/o| = (/|. This partition 
is then used to derive a Dyson type equation 



PqGqPo — 



Po 



PoG^Po - PoG^QodQoG^Po 



(6) 



with Go = G(lu) and G\ = (QqG 1 Qq) ; from this 
equation one readily obtains 
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where |/0 = Q U\f ) and (A| = (f \UQ . The pro- 
cedure can now be repeated for the diagonal element 
(/i I G\ | /i) appearing in Eq. (Q); from successive 
iterations a hierarchy of left (/„| and right |/„) states 

if )(T 1 

and corresponding projection operators P, — ' A 
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constructed according to the recursion scheme 

\fn) = U\f n -i) - a n _i|/ n _r) - 6^_x|/ n -2), 

(fn\ = (fn-i\U - a„_i(/ n _i| - bl^iU-zl (8) 



with (f n \ = \f n ) = for negative n, and where 



(fn\ U\U) . 
\fn I fn) 



(fn | fn 



(fn-1 I /n-l) 



6o = l. (9) 



From these results one can easily prove that the states 
(f n \ and |/„) form a biorthogonal set, i.e (/„ | f m ) = for 
n/m and that the values of the elements a n and b n can 
all be obtained from those of the autocorrelation function 
C(n). The operators P n and Q n = 1 — P n are not self- 
adjoint in general. However, the previous analysis only 
requires them to satisfy P% — P n and Q 2 n — Q n . 

The recursive procedure leads to the following expres- 
sion for the diagonal matrix element of the resolvent as 
a continued J-fraction 



R{z) = G ff {u) 



2^2 



bfz 



b\z 2 



1 — CLqZ — 1 — (l\Z 



1 — a 2 z 



(10) 



where z = e luJ . 

From equations (||), ^ we find readily the expression 
relating this diagonal element R{z) to the corresponding 
autocorrelation function G (n) 



(11) 



n=0 



The correspondence that we have just found between the 
power series in Eq. ( |Tl] ) and the continued fraction in 
Eq. ( |lCt ) is well known in the theory of continued frac- 
tions where a theorem establishes that under general 
conditions, such a correspondence is indeed one-to-one. 
This implies that the rational function defined by the p^ 1 
approximant 
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(12) 



to the continued fraction (^o|), which is obtained by tak- 
ing a n = b n = for n > p, has the autocorrelation values 
G(0), G(l), ... C(2p- 1) as its first 2p Taylor coefficients 
at z = 0. As a matter of fact, these are the only values 
of the autocorrelation function required in the determi- 
nation of the elements a n and b n , which are needed to 
obtain from them the p^h approximant R^(z), as can 
be deduced from their definition in Eq. (||) and from the 
recursive relations in Eq. (||). 

Let us assume for a moment that the autocorrelation 
function has the following form as a sum of p exponentials 



C(n)=X), 



(13) 



In this case the continued fraction expression for R(z) 

truncates exactly at the p^ approximant; in other words 
R(z) has an exact rational representation in the form of 
Eq. (|l2]). Then, the poles of this rational function are 
the p complex numbers z" 1 , their corresponding residues 
being -c t z~ . The analytic theory of rational functions 
shows that the converse is also true [ pq j; namely, if the 
p^h approximant R^ p \z) can be computed from the first 
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2p values of an autocorrelation function C^nVthen each 



one of these 2p values satisfy exactly Eq. fll3| ) with z% 
and Cj given, as before, in terms of the poles and the 
residues of R tyP \z). By fixing the values of the 2p pa- 
rameters (zi, Cj) from the first 2p values of C(n), we are 
indeed performing an interpolation of C (n) by a sum of 
p exponentials Jl^, [l6) . 

The memory function scheme just presented gives a 
physical support to the use of Pade approximants in the 
representation of the power series for R(z) in Eq. (11). 
For instance, from the preceding analysis we have 



R( p) (z) = 



a + ■ ■ ■ + a p -\z p 1 
1 + Piz + ■ ■ ■ + P p zP 

2p-l 

C(n)z n + 0{z 2p ). 

71=0 



(14) 



Therefore R^\z) gives the [p — l,p] Pade approximant 
to R(z). In the location of PR resonances, Pade approxi- 
mants were first used by Isola in the Henon map (11| and 
later by Baladi et al. jl2| in an intermittent system. 

Unlike the method of interpolating exponentials or its 
equivalent form as a Pade approximant to R(z), which 
can only provide the values of the leading resonances par- 
ticipating in a given observable, the memory function ap- 
proach makes also possible, in principle, the calculation 
of the corresponding generalized eigenstates. Indeed, it 
is not hard to prove (this is the essence of the Lanczos 
method jll|) that in the biorthonormal basis set given 
by the states 



|$n) = (7n-l|/n-l)- 1/2 |/„-l), 
($„| = (/„_ 1 |/„_ 1 )- 1 / 2 (/„_ 1 |, 



(15) 



the evolution operator U admits the following tridiagonal 
complex-symmetric matrix representation 



/ a h 
h ai b 2 
b 2 a 2 



U 



\ 



'•• b n 



\ 



(16) 



Again, if C(n) has the form given in Eq. (|l3|) as a sum 
of p exponentials, the above matrix truncates exactly into 
a p x p block, i.e. a n = b n = for n> p; from the eigen- 
values, Zi we obtain the location of p resonances, which 
coincide exactly with that found from the poles of the 
p^h approximant R^ p \z) in Eq. (|l2]). The correspond- 
ing right (left) eigenvectors, which are written as linear 
combinations of the right (left) states of the biorthonor- 
mal basis set given before, provide a representation of 
the right (left) generalized eigenvectors associated to the 
resonances found. 

In a real dynamical system the correlation function is 
expected to be dominated by a few leading resonances, 



but many others may intervene with small contributions. 
Besides, statistical fluctuations are always present in any 
estimate of C(n) from the system orbits. Under such cir- 
cumstances truncation of the memory function recursive 
scheme is required, but one has to choose carefully the 
order p to obtain the location of these resonances with 
enough accuracy. We will come back to this important 
issue later. 



III. FILTER DIAGONALIZATION SCHEME 

Filter diagonalization is a particular procedure to solve 
the problem of fitting a signal C(n) to a sum of complex 
exponentials as in Eq. ([13]). In this approach, this har- 
monic inversion problem is reformulated as an eigenvalue 
problem for an effective evolution operator U. The orig- 
inal formulation of this scheme Jl4|] is simplified if, as in 
our case, the signal C(n) is sampled on an equidistant 
grid Jll|. Then using our notation, the matrix represen- 
tation for U is realized in the basis set of Fourier-type 
states 



I*, 



U m \f); n = l,2, 



(17) 



where the q real phases <p n belong to an interval (p m in < 
< Vmax contained in 2ir > ip n > 0. Then the gener- 
alized eigenvalue problem to solve is 



uv = svz, 



(18) 



where the U and S are complex symmetric matrices, 
whose elements are given by 



S TO „ = (* m , *„), 
U m „ = (* m , U *„) 



(19) 



In these equations the symbol (•,•) defines a complex 
symmetric inner product, i.e without complex conjuga- 
tion. The columns of the matrix V give the eigenvectors 



(20) 



and Z is the diagonal matrix with the complex eigen- 
values Zi. Both S and U usually have rapidly decaying 
off-diagonal elements, which can all be calculate from 
the correlation values C(n). Then the eigenvalues Zi 
lying near a segment (e l¥,mi ™, e IVmax ) of the unit circle 
in the complex plane may be obtained by solving Eq. 
(|l8|) in the basis set of Eq. (|l7j ) restricted to q differ- 
ent values of the real phases tp n belonging to the interval 
(fimin < < (fimax ■ The number q should be large 
enough to extract all the leading eigenvalues in that re- 
gion of the complex plane. The amplitudes Cj are then 
calculated from the eigenvectors as = (/, Xi) |P|- 

Let us show now the close relationship between mem- 
ory function and filter diagonalization schemes. On one 
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hand, note that the use of the complex symmetric inner 
product (•, •) in the latter may be avoided if one defines, 
as in the memory function approach, a left functional 
space expanded here by the states 



p-i 



(* n | = Y,e mV "(f\U m ; n=l,2,. 



(21) 



If the observable / is a real function of the phase space 
variables (the case of complex / will be discussed later) , 
then S and U matrix elements have the usual form 



U mn = (f m I U I <&„). 



(22) 



On the other hand, we have shown in the previous 
section that the memory function scheme can be cast 
as an eigenvalue problem of a tridiagonal matrix; this 
matrix is a representation of the evolution operator U 
in a biorthonormal basis set whose states are obtained 
through a recursive procedure. From equations (g) and 
( |l5| ) we note that for any order p these states can be 
written as linear superpositions 



|$ n ) = lnmU m - X \f), n = l,2,...,p, 



p 

E 

m— 1 



n = l,2,...,p, (23) 



with coefficients 7„ m obtained recursively. Therefore, 
if we abandon the recursive scheme and instead of the 
states ($„| and |$„) we take directly 

\f(n)) = U n \f), n = 0,l,...,p-l, 

(f(n)\ = (f\U n , n = 0,l,...,p-l, (24) 

to expand the left an right functional spaces we arrive 
to a generalized eigenvalue problem equal in form to Eq. 
( |l8|) and where the p x p XJ and S complex symmetric 
matrices are now 



S mn = (f\U 



m—1 t jn—1 



l \f) = C(m + n-2), 



U mn = (flU^UU^lf) =C(m + n- 1). (25) 

The columns and rows of the matrix V give respectively 
the right and left eigenvectors 



\Xn) = Y, V mn \f(m-1)), 

(Xr, 



m—1 
V 

E 

m—1 



^ V m „(/(m-l)|. 



(26) 



The variational approach followed by Blum and Agam 
to locate the leading PR resonances in the standard 
and perturbed cat maps is a 4 x 4 implementation of this 
eigenvalue problem. 



Notice that the two basis sets defined in equations (|21 
and (|24| ) are particular choices of the general form given 
in Eq. d23|), i.e j m n = <W, and j mn = e lmiPn , respec- 



tively. Therefore, both eigenvalue problems can be made 
equivalent if there exists an invertible linear transforma- 
tions between states |<J>„) and \^n)- This occurs if q, the 
number of ip n phases in Eq. (|l7|), is chosen q = p, and 
the values e ltpn sample conveniently the whole unit circle. 
The number of correlation points C(n) required to set up 
the corresponding pxp eigenvalue problems is always 2p. 

The methods described in this and previous sections 
can be easily generalized to cross-correlation signals 
Cgf( n ) = (^(O) I f(n)). In this case, the states expanding 
the right and left eigenstates are respectively 



|/(n)> = EP|/(0)>, 
(g(n)\ = (g(0)\U n , 



(27) 



which lead to the same kind of generalized eigenvalue 
problem. As a matter of fact, the filter diagonalization 
analysis of a complex observable / corresponds in our no- 
tation to a cross-correlation function C g f(n) with g = f*. 
If a full multichannel cross-correlation matrix C'ij(n) is 
available for a set of observables /, the basis set expand- 
ing the right and left spaces may be obtained as direct 
sums of those corresponding to each observable fi (see 
also @, [§). 

With the results in this section we conclude in our 
goal of finding a unified physical framework for certain 
schemes followed to find PR resonances which make use 
of techniques such as Pade approximants and diagonal- 
ization of small dynamically adapted eigenvalue prob- 
lems. This framework clearly allows for a better analy- 
sis of issues such as convergence and numerical stability, 
which are relevant to establish the reliability of this nu- 
merical approach in locating PR resonances. These and 
other related aspects will be discussed in the following 
sections. 



IV. INTERPOLATION BY EXPONENTIALS AS 
AN EIGENVALUE PROBLEM 



In Section II we have already anticipated that when 
we truncate the hierarchy of equations appearing in the 
memory function scheme we are in fact carrying out an 
interpolation of the correlation function by a finite sum 
of exponentials jig ]. In this section we will proceed to 
give a general proof of such statement, namely we will 
show that the method of interpolating exponentials can 
be recast as a generalized eigenvalue problem equivalent 
to that of the preceding section. 

Suppose we want to solve the harmonic inversion prob- 
lem by interpolating a given set of correlation data C(n) 
by the sum of p complex exponentials J2i=i c i z ?> where 
the unknown parameters are the p amplitudes Cj and the 
p numbers Zi . If we use the 2p first values of C (n) , those 
parameters can be obtained, in principle, from the solu- 
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tion of the set on non-linear equation 

p 



C(n) = ^ Cl z™, n = 0,l,...,2p- 1. (28) 



Using these 2p correlation values we can obtain the pxp 
matrices S and U defined in Eq. (^H|) . From this equation 
and the set ( |28] ) we derive the following expressions to be 
satisfied by these matrices 



S = W W, 

u = wzw, 



(29) 



where W, and Z have matrix elements, Wy = c}J 2 z{ 1 , 
Zy = ZiSij, and W T denotes the transpose matrix. We 
finally obtain from Eq. (p9[). 



UW 



SW _1 Z, 



(30) 



which is identical to the eigenvalue problem of the pre- 
ceding section. Since W and Z on one hand, and S and U 
on the other depend on the same number of parameters 
(2p), the eigenvalue problem in Eq. ( |30|) is equivalent to 
the system in Eq. (p8|). 

A similar eigenvalue problem can be set up such that 
it is equivalent to the interpolation of the elements of a 
multichannel cross-correlation matrix Cy(n) = (fi(0) | 
fj(n)) by a sum of the same p exponentials, i.e. 



C ij( n ) = ^2 c ij,kZk, 71 = 0, 1, ... , 2pij - 1, 



k=l 



i,j = 1, 2, ...,q, 



(31) 



with the restriction cy^ = aikdjk- 

The existence of solution to the problem in Eq. j30| ) 
may be guaranteed by requiring that both S, and U be 
nonsingular matrices. A singular S matrix would reveal 
linear dependences in the basis set, and this will always 
occur when the order p, i.e. half the number of correla- 
tion data used, is larger than the true number of expo- 
nentials required to expand C(n). A singular U would 
imply a zero eigenvalue Zj, which is usually linked to a 
singular S matrix. One can get rid off these linear depen- 
dences by reducing the number of correlation data and 
thus the number of states until S becomes nonsingular. 
In the recursive method of the memory function scheme, 
if the (p + l)x(p+l)S matrix becomes singular for 
p = Po then the coefficient b Po vanishes, thus truncation 
occurs naturally. The same results are obtained if one 
performs a singular value decomposition of S JlJ|. In 
practice, in a real dynamical system, although the cor- 
relation function may be expected to be dominated by a 
few leading resonances, many others may intervene with 
small contributions. Besides, statistical fluctuations are 
always present in any estimate of C (n) from the system 
orbits. Under such circumstances there may not exist 
an order p such that b p vanishes exactly or S becomes 
exactly singular. In very favorable situations one may 



identify an optimal order p if a sharp decrease of the 
magnitude of either b p or the (p + 1) x (p + 1) determi- 
nant |S| is observed at p = p - We will see later that 
in chaotic dynamical systems such a behavior is not gen- 
erally found, and one should proceed very carefully. In 
any case, the optimal order is usually much smaller than 
the number of correlation data which can be determined. 
This fact allows for the implementation of methods to 
reduce the statistical noise in C(n), so that the 2p val- 
ues used to set up the eigenvalue problem become more 
accurate. The standard procedure is to perform a least 
squares fit. 

A general scheme for this procedure can be designed 
to analyze a full cross-correlation matrix; however for 
simplicity, we only present here the case of a single cor- 
relation signal. Our goal is then to perform an optimal 
fit of q correlation values by a sum of p exponentials, 
with q^> p. There are different solutions to this problem 
depending on the error function £ to minimize fl21j| . Here 
we will make the choice 



9-1 



n=0 



(32) 



where C(n) is the known correlation signal. 

Let us show now that this non-linear problem can be 
again reformulated as a linear one in which an eigenvalue 
problem must be solved self-consistently. By minimizing 
£ with respect to the parameters Cj and z% we obtain the 
system of equations 



^z"e„ = 0, i=l,2,...,p, 



n=0 



q-1 



^2nc i z?- 1 e n = 0, i=l,2,...,p, (33) 



n=0 



where 



e n = C(n) - 



(34) 



is the error between the known correlation values and 
their fit. We now rewrite the system ( |33] ) in the form 

2p-l 

^ %i £n tLij i l;2,...,p, 

n=0 

2p-l 

22 nc i z"~ 1 £ n = v i , i = l,2,...,p, (35) 



where Ui = X] 



n=2p z " £ ™ and v i 



y 



9-1 
n—2p 



naz. 



The solution to this system can be performed in the 
following self-consistent way. Suppose we have an ini- 
tial guess c| ^ and zj°^for the 2p parameter q and z*, 
obtained, for instance, from the solution of the pxp 
eigenvalue problem set up with the first 2p values of the 
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known correlation function C(n). We use this guess to 
calculate 



»=i 



n = 0, (36) 



and the terms Ui and vi and solve system (35) to find the 
first 2p errors e„ = e„\ With these values we obtain a 
corrected form of the correlation function 



C« +1) = £#> + V [C(n) + sip - C®] , 

n = 0,l,...,q-l, (37) 

where j = in this initial step, and r\ is a parameter 
conveniently chosen. With the first 2p values of C„ we 
reset and solve the p x p eigenvalue problem from which 
we extract an improved guess cj 1 ^ and zf^ ; the procedure 
is now iterated until the fix point of the map defined by 
Eq. ( |37| ) is reached. This will require a convenient choice 
of the parameter 77. We obtain in this way the best fit in 
the sense of the error function in Eq. ( |32"| ) and from the 
solution of the last eigenvalue problem the best eigenval- 
ues and eigenvectors. One may start with a small value of 
the number of exponentials p and increase systematically 
this number fulfilling always the condition p <C q; an op- 
timal p = p would be the maximum order for which the 
fixed point of (|3?1 ) can be numerically reached. 

Let us finish this Section with a brief discussion on the 
physical interpretations of the left and right eigenfunc- 
tions derived numerically from the schemes proposed in 
this work. Let us suppose first that the observable cho- 
sen / has a finite expansion in one or in both of the basis 
sets corresponding to the left or right generalized eigen- 
vectors of the FP operator. Then it is straightforward 
to see that an exact truncation order p = p e exists, p e 
being the number of terms in the shortest of the two ex- 
pansions. Therefore, only the p e generalized eigenvectors 
participating in this shortest expansion can be obtained 
form these numerical schemes. More accurate correla- 
tion values will provide better approximations to the ex- 
act eigenvectors. The observable / being an ordinary 
function, this situation can occur only if the eigenvectors 
of the finite expansion are themselves ordinary functions 
and not distributions. We will present an example of this 
case in the following Section. 

However, in a general chaotic dynamical systems, both 
the left and right generalized eigenvectors of the FP op- 
erator will be complicated distributions H of the phase 
space variables. In this case the two expansions of the 
observable function / will not be finite and therefore an 
exact truncation order p e will not exist. As already dis- 
cussed, we can still obtain a finite optimal order p Q , which 
will be determined by the accuracy of the correlation val- 
ues and by the numerical precision of the computations. 
The necessary choice of a finite p imposes a time cutoff 
to the evolution of the initial state and thus a limit to the 
phase-space resolution with which the generalized eigen- 
vectors can be obtained. Therefore, this class of methods 



provides in this case smooth representations of the exact 
eigendistributions. This will be also illustrated in the 
next Section. 



V. NUMERICAL EXAMPLES 

In this section we will illustrate the use of this new 
class of methods by locating the leading PR resonances of 
some simple chaotic maps. We first choose the Bernoulli 
map 



(mod. 1), 



(38) 



for which the resonance values Zi and the generalized left 
Xi and right Xi eigenstates of the FP operator are known 
analytically. These are, respectively (see Ref. B) 



1 



0,1,2, 



(39) 



and 



Xo = 0(x)9(l-x), 

Xi = (-) i - 1 [s^Hx-^-S^ix) 
_ Bj(x) 

Xi I 1 



(40) 



where S l (x) represents the fi 1 derivative of the Dirac 
delta distribution and Bi(x) are the Bernoulli polyno- 
mials. For the observable f(x) we take the polynomial 
function 



(41) 



and the corresponding normalized autocorrelation func- 

f 1 dxf(x)U n f(x) 

tion Cm) = r i can be computed analyti- 

J dxf 2 (x) 

cally from the above spectral decomposition giving 



C(n) 



14 
15' 



45 



45 



(42) 



Therefore the three resonances Zi, i = 1,2,3 [Eq. (p9[)] 
contribute to this correlation function. Of course, in a 
general practical situation we will not know the exact 
correlation function but only an estimate of it, which is 
usually obtained from the system orbits; in such case 
some statistical error will be always present in C (n) . We 
simulate these statistical fluctuations by adding to the 
exact C(n) in Eq. ( pi] ) a random noise e with Gaussian 
probability distribution 



P(e) 



1 



27TCT 



exp 



1 



Thus 



C r (n) =C(n)+e(n), 



(43) 



(44) 
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and proceed to illustrate the effect of the noise standard 
deviation a in the implementation of the methods that 
we have proposed in this work to locate PR resonances. 

A first parameter to determine is the optimal order p , 
i.e the number of resonances contributing significantly to 
the correlation function. As mentioned in Section [V, one 
may identify an optimal order p if a sharp decrease of 
the magnitude of either the memory function coefficient 
b p [Eq. (||)] or the determinant of the (p+1) X (p+1) cor- 



relation matrix S(w + 1) [Eq. (25)] is observed at p = p a . 

In figures |l| and g we plot respectively, as a function of 
p, the magnitudes of b v and \S(p + 1)| obtained from the 
noisy correlation function C r (n) for the Bernoulli map. 
We clearly observe in Fig. [j] that at the optimal or- 
der p = p = 3 the coefficient bp is very sensitive to er; 
the same sensitivity is observed in |S(p+ 1)| for p > 3 
(Fig. ||). Besides, we observe that if the error in the 
correlation function is not small enough it may wash out 
completely the expected sharp decrease in both bp and 
|S(p + 1)| at p = p ', the determination of an optimal or- 
der would therefore require quite precise correlation val- 
ues (cr < 1CT 4 ). 

If instead one makes use of the self-consistent least- 



squares method described in Section IV, p would be the 
order giving the minimum value of the error function £ 
defined in Eq. (j32|); in general, p coincides with the max- 
imum order p for which self-consistency is achieved, i.e. 
for which the fixed point of Eq. ( |37| ) can be numerically 
reached. In the present case, for a < 10~ 2 this occurs 
always at the correct value p — 3, which indicates the 



o 



-12 




3 
P 



FIG. 1: Magnitude of the memory function coefficient bp [Eq. 
(§)] as a function of the order p. These values were obtained 
for the Bernoulli map from a noisy correlation function for the 
observable f(x) = x 3 — i; the parameter a gives the standard 
deviation of the gaussian random noise. 
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FIG. 2: Magnitude of the determinant of the correlation 
matrix | S(p + 1)| [Eq. (jiH))] as a function of the order p. These 
values were obtained for the Bernoulli map. Other details are 
as in Fig. 1. 



ability of the least squares method to reduce the negative 
effects of the correlation noise. 

Let us now illustrate the effect of a in the resonance 
locations. In Table | we present the resonances derived 
from the solution of the eigenvalue problem defined in 
Section which was set up for the noisy correlation 
function C r (n); different values for a and the order p 
were used. We deduce from these results that the lead- 
ing resonance z\ is quite robust against noise and non 



TABLE I: Resonance locations of the Bernoulli map. They 
were obtained from the solution of the eigenvalue problem in 
Section which was set up for the autocorrelation function 
of the observable f(x) = x 3 — \. The parameter a is the 
standard deviation of the gaussian randon noise added to the 
exact correlation values [see Eq. (ftij)], and p is the order of the 
method, i.e. the expected number of resonances contributing 
to the correlation function; its correct value is p = 3. 



a 


p 




Z2 


23 




2 


0.4909 


-03094 




i(r 4 


3 


0.4943 


0.0881 + 0.3520i 


0.0881 - 0.3520i 




4 


0.5370 


0.4733 


-0.1207 




2 


0.4909 


-0.2891 






3 


0.4998 


0.2347 


0.1379 




4 


0.5001 


0.2659 


0.0949 




2 


0.4909 


-0.2889 




10~ 8 


3 


0.5000 


0.2499 


0.1251 




4 


0.5000 


0.2502 


0.1251 
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TABLE II: Resonance locations of the Bernoulli map. They 
were obtained with the self-consistent least-square method 
described in Section [iy|. The optimal order p = 3 given by 
this method was used. Other details are as in Table I. 



a 


21 


22 


23 


10" 4 


0.5003 


0.2697 


0.0850 


1(T 6 


0.5000 


0.2505 


0.1243 


10 -8 


0.5000 


0.2500 


0.1250 



optimal choices of the order p. The accurate location of 
other two resonances require, however, a better selection 
of p and a more converged correlation function estimate; 
for instance, the Z2 and values obtained from the most 
noisy correlation (a — 1CP 4 ) are completely wrong. 

The resonance locations obtained from the self- 
consistent least-squares method are given in Table |J. 
As we already know, this method provides as the op- 
timal order p the correct value p — 3. The method 
gives also significantly better resonance eigenvalues even 
in the case of large correlation noise (a = 10~ 4 ), where 
the previous approach failed completely. In conclusion, 
when statistical errors are present in the estimated cor- 
relation function, of all the methods discussed in this 
work, the least-squares scheme provides the most accu- 
rate resonance values by reducing the negative effects of 
the correlation errors. In general, two factors increase 
the reliability of a resonance location obtained in this 
way: a larger stability and a larger contribution in the 
observable chosen. 

The resonance generalized eigenvectors can be also de- 
termined from these schemes. For the Bernoulli map 
the calculated right eigenvectors are approximations to 
the Bernoulli polynomials participating in the observable 
chosen; a better estimate of the eigenvalue gives a bet- 
ter approximation to the corresponding Bernoulli poly- 
nomial. However the left eigenvectors provided by these 
schemes are wrong representation of the true ones. The 
reason for such an incorrect result may lie in the observ- 
able chosen and in the very different nature of the exact 
left and right generalized eigenvectors: while the right 
ones are functions, the left ones are distributions; this 
asymmetry is a consequence of the non-invertibility of the 
map. Therefore, the expansion of our observable (a poly- 
nomial) in the right eigenvector (Bernoulli polynomial) 
has a finite number of terms while this number is infinite 
if such expansion is performed in the left ones. When 
these two numbers are different, as in our case, we know 
from the final discussion in Section |^ that the present 
schemes will provide a representation for the states and 
resonances of the shorter expansion. 

As a second example we will consider the standard map 
@ 

— x n t" y n , 
K 

Vn+i = Vn + 7T" sin {2nx n+1 ) , (mod 1) (45) 

Z7T 



This map is area-preserving and invertible. While for 
K = it is integrable, for K ^ presents chaotic phase 
space regions which become more extensive as K in- 
creases. The system follows the route to chaos known 
as overlapping resonances |p2| , |23| . We take the value 
K = 10 in our calculations. 

The exact knowledge of the PR resonances is not gen- 
erally possible for this system. Blum and Agam pro- 
posed a variational approach to locate the four leading 
ones, which as we are already aware, is &p = 4 implemen- 
tation of the eigenvalue problem derived in this work. For 
comparison, we will take the observable chosen by these 
authors, i.e. 



f(x,y) = exp(i2? 



(46) 



whose autocorrelation function may be decomposed as 
the sum of the autocorrelation functions for cos 2irx (C c ) 
and sin27rx (C s ) 



C(n)=C c (n)+C s (n). 



(47) 



This decomposition is possible because the standard map 
has an inversion symmetry with respect to the point 
(x = |, y = |); therefore, the generalized eigenvectors 
are either symmetric or antisymmetric under this trans- 
formation. The symmetric ones will participate in C c (n), 
while the antisymmetric ones will do in C s (n). 

Both C c (n) and C s {n) have been determined here from 
the system orbits with an uncertainty a < 10~ 5 . We have 
proceeded next to implement the self-consistent least- 
squares method of Section IV to find the leading PR 
resonances. As a first result we obtain for the optimal 
order the values p — 9 from C c (n), and p = 8 from 
C s (n). The location of the leading resonances is given 
in Table [II. Since the two correlation functions do not 



have common eigenvalues, the total optimal order corre- 
sponding to the initial observable in Eq. ( fj6| ) is p = 17; 
this number is much larger than the value p = 4 used by 
Blum and Agam. If we take this smaller order to solve the 
general eigenvalue problem, the locations that we obtain 
for the the four resonances coincide exactly with those 
found by these authors; they are also included in Table 
III. However, if we increase p just by one these values 



TABLE III: Locations of leading resonances in the standard 
map for K = 10. The first two lines correspond respectively 
to C c (n) and C s (n) autocorrelation functions. The resonances 
were obtained from these correlation functions using the self- 
consistent least-squares scheme. The bottom line gives the 
four values obtained by Blum and Agam from their p = 4 
variational approach 



/ 




Zi 




cos(27ra;) 


0.672 


-0.030 ± 0.702i 


0.332 ±0.503i 


sin(27ra;) 


-0.715 


0.150 ±0.592i 


-0.119 ±0.553i 


e* 2 ™ 


0.515 


-0.494 


-0.003 ± 0.505i 
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0.713E3 1 



0.7051 




0.638 




0.713 



0.630 

0.550 X 0.558 0.251 X 0.259 

S 0.638 





705— B^^IIIHii »m\ 0.630 

0.550 X 0.558 0.251 X 0.259 

PIG. 3: Structure of the right (right panels) and left (left 
panels) generalized eigenvectors corresponding to the real res- 
onances z% = 0.672 (top) and Zi = —0.715 (bottom) in the 
standard map (K — 10). The eigenvectors in the top panels 
are symmetric under the inversion transformation and partici- 
pate in the cosine correlation function C c ; those in the bottom 
panels are antisymmetric and participate in the sine correla- 
tion function C s - To better illustrate the intricate structure 
we have selected small regions of the available phase space. 
Absolute magnitudes are represented, with higher values cor- 
responding to brighter regions. 



change significantly, which means that they are not prop- 
erly converged. 

The corresponding left and right generalized eigenvec- 
tors are expected to be complicated distributions. In this 
case, as we have discussed in Section [lV|, the eigenfunc- 
tions derived from our numerical approach are smooth 
representations of them. In Fig. |we present the numer- 
ical left and right eigenfunctions associated with each of 
the two leading resonances of the standard map; the first 
pair is symmetric and the second one antisymmetric un- 
der the inversion transformation. To better illustrate the 
intricate structure we have selected small regions of the 
available phase space. 



VI. SUMMARY 



tions 1 11, or the diagonalization of small dynamically 
adapted eigenvalue problems [Q. In an attempt to pro- 
vide a theoretical support to these methods, we have ana- 
lyzed their connection with other numerical schemes used 
in different contexts. A first category of such schemes in- 
cludes the memory function techniques used in the 
general theory of relaxation, and the methods related to 
this approach such as those based on the use of contin- 
ued fractions and Pade approximants. In a second cat- 
egory we have also considered the methods of the filter 
diagonalization approach |14, [l5), which is a particular 
formulation of the harmonic inversion of a time signal as 
an eigenvalue problem. 

The analysis of these schemes led us to a theoretical 
framework in which all them become equivalent formu- 
lations of the same problem: the location of the leading 
resonances contributing to a give time correlation func- 
tion. The most convenient of these formulations is as an 
eigenvalue problem from which one can obtain not only 
the poles, i.e. the resonance locations, but also a smooth 
representation of the generalized eigenvectors. Besides, 
we have proved that all these procedures are also partic- 
ular linear formulations of the non-linear numerical prob- 
lem known as interpolation by exponentials |jl2| , This 
connection has provided us with new tools to perform a 
better analysis of issues such as convergence and numer- 
ical stability. From such analysis improved schemes may 
be designed, as the one that we have proposed based on 
a least-squares fit. 

We have illustrated the use of this class of methods in 
two chaotic maps: the Bernoulli map and the standard 
map. In the first example the nature of the right eigen- 
vectors, which are known to be the Bernoulli polynomi- 
als, allows for finite expansion of certain observables like 
the one chosen for illustration. Then an accurate enough 
correlation function provides good approximations to the 
PR resonances and their corresponding right eigenvec- 
tors. The least squares method improves significantly 
these results in the case of less accurate correlation val- 
ues. 

In the standard map, for which the exact knowledge 
of the PR resonances is not generally possible, the more 
elaborated least squares method provides resonances val- 
ues which appear to be significantly more converged than 
previous results 0. Nice smooth representations of the 
eigendistributions are also obtained. These indicate the 
relevant phase space regions involved in the dynamics. 



Among the schemes followed in the literature to locate 
the leading Pollicott-Ruelle resonances of the Perron- 
Frobenius operator in chaotic dynamical systems, there 
are a few examples of remarkable simplicity, which in- 
volve either the use use of Pade approximants to perform 
the analytical continuation of the spectral density func- 
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